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In this paper we apply digital in-line holography to image opaque objects 
through a thick plano-concave pipe. Opaque fibers and opaque particles are 
considered. Analytical expression of the intensity distribution in the CCD sensor plane 
is derived using generalized Fresnel transform. The proposed model has the ability 
to deal with various pipe shape and thickness and compensates for the lack of 
versatility of classical DIH models. Holograms obtained with a 12 mm thick 
plano-concave pipe are then reconstructed using fractional Fourier transform 
(FRFT). This method allows us to get rid of astigmatism. Numerical and experimental 
results are presented. © 2012 Optical Society of America 

OCIS codes: 090.0090, 070.0070, 100.0100 
1. Introduction 

Digital in-line holography (DIH) is a recognized optical technique for flow measurements in trans- 
parent liquid media. This method is widely used in various domains such as fluid mechanics [TJH] 
where the flow is seeded with small particles, or in the biological microscopic imaging field J3JH]. 
Nevertheless, a transparent pipe can be viewed as an optical system that can introduce 
aberrations such as astigmatism [5]. These aberrations make difficult to image objects 



in the pipe. Compensation for these unwanted effects has been widely investigated. 
De Nicola et al. |6j proposed a numerical method to compensate for anamorphism. By 
modifying the chirp function and the spatial frequency term of the Fresnel integral 
the authors managed to compensate for a severe anamorphism brought by a reflex- 
ion diffraction grating. Numerical methods for astigmatism compensation have been 
proposed in Refs. [6H8]. Here the authors used a modified chirp function where two 
different propagation distances are considered. Astigmatism can also be compensated 
for by optical means. For instance, the authors of [9j used index matching to mini- 
mize astigmatism: the thin-lens like effect of an ampule filled with contaminants was 
investigated to determine index matching parameters for astigmatism compensation. 
However, the pipe used to carry a flow must be considered as a thick and cylindrical op- 
tical system. Recently, an analytical solution of the scalar diffraction produced by an opaque disk 
illuminated by an elliptical, astigmatic and Gaussian beam under Fresnel approximation has been 
proposed [10] . In this publication, the astigmatism was introduced and controlled by using a thin 
plano-convex cylindrical lens. Using fractional Fourier transformation (FRFT), authors managed 
to retrieve a correct image of an opaque disk from astigmatic hologram. It was shown that FRFT 
is therefore well suited for these studies where classical methods such as Fresnel transform [11] or 
Wavelet transform [T^] fail. 

The astigmatism introduced by thin lenses is thus well controlled. However, these DIH models can 
not be applied to hugger pipes. As a matter of fact, the optical thickness of these has to be taken 
into account. 

In this paper, the aim is to apply DIH to thick pipe systems. In the first part of the paper, the 
intensity distribution, in the CCD plane, of the diffracted field produced by an opaque object located 
in a pipe is calculated. A method of Gaussian functions superposition is used to describe 
the object function. In the second part, definition of the FRFT is recalled and we demonstrate 
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its ability to reconstruct holograms recorded in such systems. Finally simulations and experiments 
are performed to illustrate our results. 

2. In-line Holography through pipes 

Holography aims to record, on a CCD camera without objective lens, the intensity distribution of the 
diffraction pattern of an object illuminated by a monochromatic continuous wave [13] . The numerical 
and experimental set-up is represented in Fig. ([TJ. The incident Gaussian beam propagates in free 
space over a distance z p and illuminates the pipe. Here, the pipe is modeled as two thick lenses. Their 
thickness is denoted by e. The opaque object is located between these two lenses: at a distance 
S from the first thick lens, and at a distance Zj from the second lens. The CCD sensor 
is located at the distance z next the pipe and records the intensity distribution of the diffraction 
pattern. 

2. A. Intensity distribution in the CCD sensor plane 

In this part, we consider the propagation of a Gaussian beam through our optical system. In the 
beam waist plane, located at a distance — z p from the pipe, the Gaussian beam, denoted G, is defined 



where uj is the waist width and (/i, v) are the coordinates in the beam waist plane. 

The propagation of the Gaussian beam through the pipe, to the CCD sensor is decomposed into 
two steps. The first step is the propagation of the illuminating beam (i.e. G(n,v)) from the beam 
waist plane to the front of the opaque object. Using the generalized Huygens-Fresnel integral [HMl6| . 
the amplitude of the field in the pipe, denoted G\, is defined by: 



by: 
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A\ ,v , B^ v , D\' y are given by M^' v matrices of Eq. ([55)1 in App. (|A"|). To obtain the relation of 
the intensity distribution, the linear canonical transformation is used. For the basic definitions and 
properties of the linear canonical transformation, we refer to |17j . This transformation is a parame- 
terized integral operator of parameters A, B, C and D. Each parameter represents a coefficient of a 
transfer matrix denoted M (Cf appendix |A"|. In our case the Fresnel transform is considered [T51ll9j . 
It should be noted that the values of the coefficients arc different in both £ and rj direction. This is 
due to the cylindrical geometry of the pipe. The distance E\ = z p + e + 5 corresponds to the waist- 
object distance. Using the same formalism, the amplitude in the CCD plane, denoted G 2 (x,y), can 
be written as: 
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where the parameters A 2 ,v , B 2 ' y , D 2 ' y are given by M 2 ' v matrices of Eq. (|39|) . The distance E 2 = 
Zi + e + z is the distance between the object and the CCD sensor. 

The spatial transmittancc of the opaque 2D-object is defined by [I — T(£, 77)]. Here T(£, 77) can 
be expressed as a superposition of Gaussian functions [2U1I2T] . such as: 
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This expression permits to deal with non symmetrical optical systems and elliptical opaque ob- 
jects p2H25] . This Gaussian decomposition is very convenient; as a matter of fact it allows us to 



establish analytical expression of the diffracted pattern in the CCD plane. It also allows us to 
simulate holograms of 3D opaque objects whose 2D projection is elliptic or circular |26j 
(e.g. spheroids, fibers ...). 

The object parameters are illustrated on Fig. ([2]). The coefficients a and b are the object radii 
within 77 and £ axis respectively, 9 is the angle between elliptical aperture principal axis and £ axis. 
The Ak and B/~ coefficients are determined by numerical resolution of the Kirchhoff equation |21) . 
Let R e u = b/a representing the particle cllipticity. Considering R e u = I lead to the simulation of a 
circular particle, whereas R e u ^ 1 is used to simulate elliptical particles. The particular case of the 
opaque fiber, parallel to the pipe axis, is obtained when R e u — > 0. 

In further developments, 9 will be considered to be tt/2. With this assumption T (£,77) becomes: 



To simulate our opaque objects, N is fixed to 10. From Eq. ©, G2 (x, y) is split into two integrals 
denoted R (x, y) for the reference beam and O (x, y) for the object beam so that : 
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The functions R (x, y) and O (x,y) are respectively given by: 
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Values of the different parameters of Eqs. (fTTj) and (|T2|) arc defined in App. (|Bl) . 

After theoretical developments the intensity distribution, denoted I(x,y), recorded by the CCD 



sensor is: 



I (x, y) = G 2 (x, y) G 2 (x, y) = ^ (\R\ 2 - 2K {RO} + \0\ 2 ) , 



(13) 



where the upper bar denotes the complex conjugate and 3? the real part. The square modulus |i?| 2 
corresponds to the directly transmitted beam whereas \0\ 2 is associated with the diffracted part of 
the beam. 

To illustrate Eq. (|T^| , the optical set-up given on Fig. (JTJ) is considered. We investigate the diffrac- 
tion pattern obtained with an opaque fiber parallel to the pipe axis. From the Gaussian decom- 
position of the object function, opaque fibers, parallel to the pipe axis, can be obtained taking 
Reii 0. The value of the opaque fiber diameter is 2b = 5f.8 m and a — > oo. The glass-made 
(m = f .5) pipe is filled with water of refractive index n 2 = 1.33. The beam propagates in free space 
over z p = 325 mm. The dimensions S and z are fixed to 18 mm and 23 mm respectively. 

The image of Fig. §S§ illustrates the intensity distribution of the diffracted field. The simulation 
is carried out by calculating the gray level of each pixel using Eq. ([TU)) . The size of the hologram is 
768 x 576 pixels, and pixel pitch is 11 m. To validate our simulation, an experiment, which consists in 
placing an opaque fiber in the pipe, was performed. The hologram of Fig |4]) represents the intensity 
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of the diffraction recorded with the CCD sensor using the same parameters than previously. This 
illustration reveals a good accordance between numerical and experimental diffraction patterns. To 
confirm this point, the transverse intensity profiles obtained from Fig. ([3]) and Fig. (jj) arc presented 
on Fig (|5|). Here the normalized intensity I* is plotted against the x-axis. The intensity profiles are 
calculated by cumulating the gray levels along the 77-axis (direction of the fiber) over 50 rows, which 
are figured out by the rectangular selection on Fig. (j4]). The diffraction pattern of Fig. 
(|4|) was shifted by 75 m along x-axis so that numerical and experimental results can be 
compared. This figure shows the good agreement between numerical and experimental data. 

2.B. Case of thin lenses 

We have proposed a theoretical model allowing us to deal with propagation of a laser beam through 
a thick pipe. In a former study, the effect of cylindrical lenses on the diffraction pattern of a particle 
has been investigated in details, leading to the expression of the intensity distribution in the CCD 
plane |10j . Using this approach we have a great opportunity to confirm our pipe model results. 

In the following developments we aim to compare our thick and cylindrical pipe model with a thin 
cylindrical lenses approach. The phase transformation due to a thin lens is jT0l|T8] : 



Here (xi,yi) are the coordinates in the lens plane and f x>y are the focal lengths of the lens in both 
directions. Using Huygens-Fresnel integral |18j the intensity distribution in the sensor plane is: 



Here, Rthin and O t hm can be expressed in a quite similar form than amplitude distributions given 
in Eqs. (HJ) and (TH|) . 

The transverse intensity distribution profile presented in Fig. (|6]) shows: 




(14) 
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remembering that e represents the glass thickness. This term is contained in matrices 
Mf ' y and M%' v (see Eqs. dMD and fl39|)). 

For this example: 2b = 51.8 m, a — > oo, z p = 325 mm, ri2 = 1, 5 = 18 mm, and z = 23 mm. The 
comparison with other parameter values leads to the same conclusion. As a result, we can consider 
that our results are consistent with those of Ref. [TU] and that our model is versatile enough to deal 
with various pipe shapes and thickness. 

In this section a numerical model allowing to treat thick optical systems has been presented. We 
now aim to reconstruct the image of the object by means of FRFT from the calculated intensity 
distribution of the diffraction pattern (Eq. (fTB")) '). In section [3J mathematical definition of the FRFT 
is recalled leading to the reconstruction of holograms recorded in pipe systems. 

3. Fractional Fourier transformation analysis of in-line holograms 

3. A. Two-dimensional Fractional Fourier transformation 

The FRFT is a generalization of the classical Fourier transform. This integral operator has numerous 
applications in signal processing fTT] . Its mathematical expression is the following [2TH21?] : the FRFT 
of order a x — (2a x )/iT and a y = (2a y )/ir (for x and y cross section respectively), with < \a x \ < ir/2 
and < \a v \ < 7r/2, of a two dimensional function I (x,y) is 
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Here p = x,y. Generally, the parameter s p is considered as a normalization constant. In our case, 
its value is defined from the experimental set-up according to [30]: s 2 = N^ x ■ S 2 . The number 



of samples is Ng x in both intensity distribution / {x, y) and fractional domain. 6 P is the sampling 
period along the two axes of the image. C{a p ), which is a function of the fractional order, insures the 
energy conservation law to be valid. Discretization of the FRFT kernel is performed using 
an orthogonal projection method proposed by Pei et al. |31j . 

3.B. Optimal fractional orders to refocus numerically over the object 

In our case, the numerical reconstruction can be considered as a numerical refocusing over the 
object [32]. To do that, the quadratic phase, denoted ip = arg(i?Oj, and contained in the term 
25i {i?0} of Eq. ([13]) must be evaluated. This term is composed of a linear chirp modulated by a 
sum of complex Gaussian functions. Information about the distance between the CCD sensor and the 
object is carried by the linear chirp (£?2 X „ ) whereas information about object size is carried by the 
modulation (M x . y ^ and N Xt y^). The optimal reconstruction of the image consists in compensating 
the quadratic phase terms of the intensity distribution [10] . 

The quadratic phase can be determined from Eqs. (fTT]) and (fT2]) 
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thus K {RO} can be written as: 

U{RO} = | RO | cos (iip). (21) 

The quadratic phase term contained in the FRFT kernel, denoted (p a , is given by 
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The image reconstruction is obtained by applying the FRFT to the intensity distribution of Eq. 
&■ 

& ax , ay [I{x,y)] oc ^ ax , ay [\R\ 2 + \0\ 2 ] - 2^ ax , ay [\RO\ cos?]. (23) 

The terms |i?| 2 and \0\ 2 contain no linear chirps, thus they do not influence the optimal fractional 
orders of reconstruction to be determined. Only the second term is useful for image reconstruction. 



By noting that 2 cos <p = exp (— ip) + exp (ip), Eq. (f2U)) becomes: 
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The best hologram reconstruction is reached when one of the quadratic phase term is brought to 
zero. Thus reconstruction is performed if 



tp a ± cp = 0. 



(25) 



The optimal fractional orders a° pt and a° pt are defined from Eqs. (|20[) . (|22|) and (|25j) and take the 
values: 



aZ pt = arctan 
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It should be noted that we have checked that if we consider M2 X to be free space propagation, 
leading to B^ y = Zi + e + z and D% ' v = 1, the expression of the optimal fractional order of Ref. [TU] 
are recovered. 

3. C. Numerical simulations 

Experimental set-up of Fig. ([T]) is used to perform simulations. An application of the model is 
illustrated on Fig. ([3]). Recall that we consider an opaque fiber (R e u — > 0) 2b = 51.8/im in width. 
The glass made pipe is filled with water (ri2 = 1.33). S and z are set to 18 mm and 23 mm 
respectively. Rx\, Ryi, Ryi — > oo and Rx2 = 18 mm are the first thick lens curvatures along x- and 
y-axis. The curvature radii of the second thick lens have the same values but opposite signs. With 
these parameters and owing to Eq. (|26[) the optimal fractional orders within this configuration are 

(27) 



a°f = 0.337, a° y pt = 0.273. 



The reconstruction of the fiber with FRFT can be seen on Fig. ([7]). As we can see on this figure, the 
reconstructed image is disturbed by background fringes. This effect is due to the in-line configuration 
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and is commonly known as the twin image effect |33j . 
3.D. Experimental results 

Since now, we have presented a method to simulate holograms in thick optical systems. Reconstruc- 
tion has been successfully performed thanks to FRFT. In order to validate theoretical developments 
and simulations, a glass pipe with curvatures Rx±, Ry±, Ry% — > oo and RX2 = 18 mm has been used 
for the experiments. Curvatures on the other side of the pipe are deduced from Rx\, Ry\, Rx2, Ryi 
by taking negative values. 

The image of Fig. Q represents the intensity of the diffraction pattern recorded with a 768 x 
576 px 2 CCD sensor with 11 \xm pitch. The object is a 26 = 51.8/im opaque fiber, S and z are 
approximately and respectively equal to 18 mm and 23 mm. Theoretical orders associated with this 
experiment are given by Eq. (|27|) . To perform FRFT on this experimental image, FRFT orders have 
been adjusted to obtain the best image of the fiber i. e the best contrast between the reconstructed 
fiber and the background. Thus doing, the accuracy on the fractional order is approximately 10 ~ 2 . 
Reconstruction is presented on Fig. (|8|), estimated optimal fractional orders are a x = 0.33 and 
dy = 0.27. These values are very close to the theoretical ones. 

This good accordance is confirmed by the curves of Fig. © . Here theoretical values of a° pt (solid 
line) and a° pt (dashed line) are plotted versus the distance between the pipe and the CCD sensor. 
Theoretical values corresponds to an air filled glass pipe, i5 is fixed to 10 mm and z varies from 
to 140 mm. Note that, when z = 0, a x . y ^ 0, this is due to the fact that, in this case, z represents 
the distance between the pipe and the CCD sensor instead of the distance between the object and 
the CCD sensor which is considered in |10j . Experiments performed with this configuration allow us 
to plot estimated values of a° pt (circles) and a opt (diamonds) on the same graph. As shown by Fig. 
©, estimated and theoretical values are closely linked. Theoretical representation hereby presented 
is therefore well adapted to these optical systems. 
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Until then, we have performed reconstruction with objects for which position is well 
known. We now apply our formalism to reconstruct hologram of latex beads. The 
intensity distribution of the diffracted field is represented on Fig. (|10p . Here, 100 m 
in diameter latex beads are used, z p = 325 mm (i.e. propagation distance between the 
source and the pipe) and z — 23 mm (i.e. distance between the pipe and the CCD 
sensor). This study gives us the opportunity to compare reconstruction using our novel 
approach with reconstruction using a thin lens approach. As far as every parameter 
of the experimental set-up is known, excepting 5, the reconstruction process is the 
following: for S G [0,36 mm] we calculate the optimal FRFT orders using Eq. (|26[). Doing 
so, we are able to refocus on the particles which are in the pipe. 

In Fig. ((6|), it is shown that when the pipe thickness e is close to zero, the two 
approaches are equivalent. Thus, to reconstruct the hologram of Fig. (|10[) with the thin 
lens approach, we only need to use the orders given in Eq. (|26l) with e = 0. 

Comparison between the two reconstruction methods can be made with Fig. (| 1 1|) . 
Here, the reconstruction is realized with the same value of 8 in both cases. We can notice 
that refocusing on the particle is impossible with the thin lens approach, whereas our 
thick lens approach allows a good reconstruction of the image of particles. 

Therefore, the thick lens formalism is well adapted to pipe flow studies: as we are 
able to refocus on particles in a pipe, metrologies of particles diameter is possible. 

4. Conclusion 

In this paper, an analytical expression of the scalar diffraction, under Fresnel approximations, in 
thick optical systems such as pipes is derived. As such optical systems reveals astigmatism, FRFT is 
a reliable tool to perform hologram reconstruction. Optimal fractional orders have been calculated 
leading to a satisfactory reconstruction of either numerical or experimental images. Numerical ex- 
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periments have been performed showing a good agreement with experimental results. Comparison 
between thin lens and thick lens approaches has been performed showing that the thin 
lens formalism is not well adapted to the study of pipe flows. 
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A. Appendix A : Transfer matrices of the optical system 



Each part of an optical system can be represented by a matrix. Using this principle, we can build a 
set of matrices corresponding to our experiments (Fig. {J}). After the beam shaping step, the beam 
propagates in free space over a distance z p . The associated M z matrix is: 



' 1 z n ^ 



M 2 „ = 



(28) 
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Then the beam is refracted at the interface between free space and glass (m = 1.5). The curvatures 
along both directions are Rx\ and Ry\. Mr x1 , Mr v1 matrices along x and y axis are: 
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After refraction, the beam propagates in glass (m=1.5). For this step: 
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Next step is the refraction of the beam at the interface between glass and the medium inside the 
pipe (refractive index n%). Curvature are given by Rx2 and Ry2 in x and y direction. 
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We are now in the pipe. To reach the object, we have to propagate over S. 



M 5 = 
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(32) 



Doing the same over Zi permit to reach the output of the pipe: 
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After a refraction at the interface (curvature Rx% and Rys) between the medium with refractive 
index 77,2 and glass: 
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a refraction at the interface (curvature Rx/± and Ry^) between glass (ni=1.5) and free space: 
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and free space propagation over z: 
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the whole ABCD system is described. 



B. Appendix B : Amplitude distributions R(x,y) and 0(x,y) 

In Appendix |A1 each part of our optical system has been represented with transfer matrices. Using 
this formalism allows us, under paraxial conditions, to deal with propagation of a Gaussian point 
source through the pipe. 

Intensity distribution of the diffraction pattern in the CCD sensor plane is determined by consid- 
ering two matrix systems: M^' v and M^- 

M*' v is composed of three steps : propagation in free space over z p , propagation through the first 
thick lens, propagation in a medium of refractive index n,2 over 6. It is characterized by two transfer 
matrices: 



( A x,y ,,,..< 
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M% ' v is also composed of three steps : propagation in a medium of refractive index n-i over Zi 



propagation through the second thick lens, propagation in free space over 8. Transfer matrices for 
this system are: 



Ml ' y = M z x ,\///" x M z 
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Thanks to Eqs. ([29)1 , (|30j) . (|31j) one can build the transfer matrices of the first thick lens: 



M x Li = M RX2 x M e x M Rxi , M y Li = M Ry2 x M e x M Ryi . 



(40) 



By using the same method for the second lens, we obtain: 



Ml 2 = M Rxi x M e x M RX3 , Ml = M Ryi x M e x M Ry3 



(41) 



B.A. Propagation through M^' v 

After analytical developments of Eq. ([2)1. the complex amplitude distribution in the object plane is: 
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y\ x K x K v exp 
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and i?i x are respectively the beam radii and the wavefront curvature in the particle plane 
Their mathematical expressions are: 
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B.B. Propagation through M. 
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Propagation to the CCD sensor plane is calculated thanks to the generalized Huygcns-Fresnel integral 
(see Eq. ©). 
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B.B.I. Amplitude distribution R(x,y) 



R is associated with the reference wave. After analytical developments of Eq. ([9]), the amplitude 
distribution R is: 
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B.B.2. Amplitude distribution O {x,y) 

O is the amplitude of the diffracted wave. We define wi x and i?i x 
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to simplify notations. It should be noted that 5R and 9 stand for real and imaginary part respectively. 
Thus O becomes: 
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D. 



, (50) 
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with 



■ >■■!) 



N x ,y e = — "4^7 r~2 > (51) 
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and 



1/2 



XBi 



(52) 
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Fig. 1. a) Schematic representation of the optical set-up (not to scale). Definition of the 
numerical and experimental parameters, b) Close-up of the pipe used in simulations and 
experiments 
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Fig. 2. Schematic representation of the object. 
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Fig. 3. Simulation of the diffraction pattern of a 51.8 \im opaque fiber, parallel to the axis 
of a glass pipe, n,2 = 1.33, A = 632.8ti?ti, z = 23mm, S = 18mm. 



25 




26 




Fig. 6. Normalized intensity distribution recorded by the CCD: comparison between thin 
lens and thick lens models, with e = mm. Simulation parameters are : 2b = 51.8 to, 
a — > oo, z p = 325 mm, ni = 1, S = 18 toto and z = 23 mm 
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Fig. 8. Reconstruction of the fiber image from the diffraction pattern of Fig. |4]) by FRFT 
with a° pt = 0.33 and of = 0.27. 
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Fig. 9. Comparison between theoretical fractional orders and optimal fractional orders esti- 
mated from the experimental holograms. 
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